module track

    use global
    use criteria, only: check
    use io, only: pos

    implicit none
    private

    public :: track_particle

    contains

    subroutine track_particle(i0, j0, k0, l0, rx0, ry0, rz0)

        integer, intent(in) :: i0, j0, k0, l0
        real, intent(in) :: rx0, ry0, rz0

        integer :: i, j, k, l
        real :: rx, ry, rz

        real :: s0, t0
        logical :: die
        integer :: n
        real :: s, se, sw, sn, ss, st, sb, s1
        real :: vol

        i = i0
        j = j0
        k = k0
        l = l0
        rx = rx0
        ry = ry0
        rz = rz0

        ! Initialize particle. This is the time counter in seconds (t0) and in
        ! scaled time (s0 = t0/volume).
        t0 = 0.0
        s0 = 0.0

        ! Main particle loop.
        do n = 1, int(maxsteps)
            ! Check special conditions to see if this particle should stop.
            ! We pass the initial and current position.
            call check(i, j, k, l, rx, ry, rz, t0, die)
            if (die) goto 10

            ! Current scaled time (s0) and interval (s1).
            vol = volume(i,j,k,1+mod(l-1,dldm))
            s0 = mod(t0, interval)/vol
            s1 = interval/vol

            ! Calculate travel time for each axis.
            call travel(bx(i,j,k,l), dx(i,j,k,l), rx, s0, sw, se)
            call travel(by(i,j,k,l), dy(i,j,k,l), ry, s0, ss, sn)
            call travel(bz(i,j,k,l), dz(i,j,k,l), rz, s0, st, sb)
            s = min(s1, sw, se, ss, sn, st, sb)

            ! Move particle inside box.
            call move(rx, s, s0, ax(i,j,k,l), bx(i,j,k,l), gx(i,j,k,l), dx(i,j,k,l))
            call move(ry, s, s0, ay(i,j,k,l), by(i,j,k,l), gy(i,j,k,l), dy(i,j,k,l))
            call move(rz, s, s0, az(i,j,k,l), bz(i,j,k,l), gz(i,j,k,l), dz(i,j,k,l))

            ! Update time counter.
            t0 = t0 + (s - s0)*vol
            if (s == s1) l = 1 + mod(l, ldm)  ! increase ``l`` if time expired

            ! Update particle position.
            if (se == s) then
                rx = 0.0
                if (cyclic_boundaries) then
                    i = 1 + mod(i, idm)
                elseif (i < idm) then
                    i = i+1
                endif
            endif
            if (sw == s) then
                rx = 1.0
                if (cyclic_boundaries) then
                    i = mod(idm + i - 1, idm)
                elseif (i > 1) then
                    i = i-1
                endif
            endif
            if ((sn == s) .AND. (j < jdm)) then
                ry = 0.0
                j = j+1
            endif
            if ((ss == s) .AND. (j > 1)) then
                ry = 1.0
                j = j-1
            endif
            if ((sb == s) .AND. (k < kdm)) then
                rz = 0.0
                k = k+1
            endif
            if ((st == s) .AND. (k > 1)) then
                rz = 1.0
                k = k-1
            endif

        enddo
        10 continue

    end subroutine track_particle
    
    subroutine travel(b, d, r0, s0, sa, sb)
        
        real, intent(in) :: b, d, r0, s0
        real, intent(out) :: sa, sb

        real :: v, fr0, fr1

        sa = 1e30
        sb = 1e30

        ! Velocity at point.
        v = -b*r0 - d

        if (b == 0.0) then
            if (d /= 0.0) then
                if (v > 0.0) then
                    sb = (r0 - 1.0)/d
                elseif (v < 0) then
                    sa = (r0 - 0.0)/d
                endif
            endif
        else
            fr0 = r0 + d/b
            if (v > 0) then
                fr1 = 1.0 + d/b
                if ((fr1/fr0) > 0) sb = s0 - (1.0/b)*log(fr1/fr0)
            elseif (v < 0) then
                fr1 = 0.0 + d/b
                if ((fr1/fr0) > 0) sa = s0 - (1.0/b)*log(fr1/fr0)
            endif
        endif

    end subroutine travel

    subroutine move(rs, s, s0, a, b, g, d)

        real, intent(inout) :: rs
        real, intent(in) :: s, s0, a, b, g, d

        if (s == s0) return

        if (b == 0) then
            rs = rs - d*(s-s0)
        else
            if ((rs + d/b) == 0.0) then
                rs = -d/b
            else
                rs = (rs + d/b) * exp(-b*(s-s0)) - d/b 
            endif
        endif

    end subroutine move

end module track
